/*
9.27
1.function matrix_calc_taskA,matrix_calc_taskB
2."N" means iteration times
10.2
1.merge the for loop in taskB
2.adding multithread for taskA and taskB.
3.move memory access variable out of "for loop".
10.12
second loop will be faster?

*/
#ifndef _CALC_H_
#define _CALC_H_

#ifdef __cplusplus
extern "C"{	
#endif


#include<omp.h>
#include<ctime>
#include<stdio.h>
typedef struct {
	int* 			rowArray;
	const int* 		rowOffset;
	int 			rowArraySize;
	const int* 		columnIndice;
	const double* 	S;
	const double*	valueNormalMatrix;
	double* 		Id;
}TaskMatrixInfoA;

typedef struct {
	const double * 	valueSpiceMatrix;//N*2
	const int* 		rowOffset; 		//n+1
	const int* 		columnIndice; 	//N

	double* 		A;	//N
	double*		 	S;	//nx1
	double*			R;	//nx1
	double* 		H;	//nx1
	const double* 	D;	//N*2
	double* 		IC;	//nx1
	double* 		IG;	//nx1
	double 		alpha;
	int * 			rowArray; //kx1
	int 			rowArraySize;
	void* 			hdl;  //? what's this?
}TaskMatrixInfoB;

#define thread_num 16
void matrix_calc_taskA(TaskMatrixInfoA** listDataList,int N)
{
	clock_t start,end;
	for (int i = 0 ;i < N; ++i)
	{	
		int row_array_size= listDataList[i]->rowArraySize;	
		TaskMatrixInfoA * list_Data_List=listDataList[i];
		const double *valueNormalMatrix = list_Data_List->valueNormalMatrix;
		const double *S = list_Data_List->S;
		const int *columnIndice = list_Data_List->columnIndice;
/*		
		const int *rowArray = list_Data_List->rowArray;
		const int *rowOffset = list_Data_List->rowOffset;	

		int temp_offset1[row_array_size];
		int temp_offset2[row_array_size];
	start = clock();

//	#pragma omp parallel for num_threads(8)schedule(guided)
//		for (unsigned int i=0;i<row_array_size;i++ ){
//			temp_offset1[i]= rowOffset[rowArray[i]];
//			temp_offset2[i]= rowOffset[rowArray[i]+1];
//		}
//	end = clock();
//	printf("%d \n",end-start);
	start = clock();
*/
//	#pragma omp parallel for num_threads(8)schedule(guided)
	int j;
	#pragma omp parallel for num_threads(2) private(j)
		for (unsigned int it = 0; it < row_array_size; ++it	)
		{
			const int node = list_Data_List->rowArray[it];
			int colum_num = list_Data_List->rowOffset[node + 1];
//			double &id=list_Data_List->Id[node];
			
			double t = 0;
			for ( j = list_Data_List->rowOffset[node];j < colum_num;++j)
//			for (int j = temp_offset1[it];j < temp_offset2[it];++j)
			{
				unsigned short colNum = columnIndice[j];
				t += valueNormalMatrix[j]*S[colNum];
			}
			list_Data_List->Id[list_Data_List->rowArray[it]] += t;
		}
//	end = clock();
//	printf("%d \n",end-start);
	}
}

//Task B interface
void matrix_calc_taskB(TaskMatrixInfoB** listDataList,int N)
{
//int coreNum = omp_get_num_procs();
omp_set_num_threads(8);
#pragma omp parallel
	{
#pragma omp for
	for (int i = 0;i<N;i++){
		TaskMatrixInfoB * list_Data_List = listDataList[i];
		int row_array_size = list_Data_List->rowArraySize;
		const double alpha = list_Data_List->alpha;
		double* A = list_Data_List->A;
		const double* S = list_Data_List->S;
		const int* columnIndice=list_Data_List->columnIndice;
		const double* valueSpiceMatrix=list_Data_List->valueSpiceMatrix;

		for(int it = 0;it < row_array_size;++it)
		{
			int row = list_Data_List->rowArray[it];
			const int k1 = row * 2;
			double current = list_Data_List->D[k1];
			double charge = list_Data_List->D[k1+1];
			int colum_num = list_Data_List->rowOffset[row + 1];
									
			double &IG = list_Data_List->IG[row];
			double &IC = list_Data_List->IC[row];
			
			for (int p = list_Data_List->rowOffset[row];p < colum_num;++p){
				//int col = list_Data_List->columnIndice[p];
				int col = columnIndice[p];
				const int k = p * 2;
				double cond =valueSpiceMatrix[k];
				double cap  =valueSpiceMatrix[k+1];					
				double ig = cond * S[col];
				double ic = cap * S[col];
				//(4)(5)
				//list_Data_List->IG[row] += ig;
				//list_Data_List->IC[row] += ic;
				IG += ig;
				IC += ic;
				//(7)(8)
				current -= ig;
				charge -= ic;
				//(9)
				//list_Data_List->A[p] = cond +  cap * listDataList[i]->alpha;
				A[p] = cond +  cap * alpha;
			}
 			//(7)(8)
			list_Data_List->R[row] = current; 
			list_Data_List->H[row] = charge;
		} 
	}
	}
}

#ifdef __cplusplus
}
#endif

#endif
